//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Quadruple_Range_2D.cc
 * \author Jeremy Roberts
 * \date   06/19/2011
 * \brief  Quadruple_Range_2D member definitions.
 * \note   Copyright (C) 2011 Jeremy Roberts
 */
//---------------------------------------------------------------------------//
// $Rev:: 141                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-09-14 19:28:23 +0000 (Wed, 14 Sep 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <cmath>
#include <cstdio>

#include "angle/Quadruple_Range_2D.hh"
#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"

namespace slabtran
{

//---------------------------------------------------------------------------//
// CONSTRUCTOR
//---------------------------------------------------------------------------//
/*!
 * \brief Constructor.
 *
 * \param sn_order   number of azimuthal angles
 * \param norm       normalization of the quadrature weights
 *
 */
Quadruple_Range_2D::Quadruple_Range_2D( int    sn_order,
		                                double norm )
: Base(sn_order, norm)
{
	using util::soft_equiv;

	Require (d_sn_order >= 2);
	Require (d_norm > 0.0);
	Require (d_sn_order % 2 == 0);

	// size angles, 2D only
	d_angles.resize(2);

	// set pointers to mu, eta, xi
	Base::Vec_Dbl *mu  = &d_angles[MU];
	Base::Vec_Dbl *eta = &d_angles[ETA];

	// number of azimuthal angles (per octant)
	int num_phi = d_sn_order;

	// number of polar angles (per octant)
	int num_theta = num_phi / 2;

	// number of octants
	int num_oct = 4;

	// total number of angles in the quadrature
	int num_ang = num_theta * num_phi * num_oct;

	// Force the direction vectors to be the correct length.
	mu->resize(num_ang);
	eta->resize(num_ang);
	d_wt.resize(num_ang);

	// compute mu and eta for the first quadrant
	double sintheta = 0;
	double wtheta   = 0;
	double sinphi   = 0;
	double cosphi   = 0;
	double wphi     = 0;
	int k = 0;
	for ( int i = 0; i < num_phi; i++ )
	{
		cosphi = get_phi( num_phi, 0, i );
		sinphi = std::sqrt( 1 - cosphi*cosphi );
		wphi   = get_phi( num_phi, 1, i );
	    for ( int j = 0; j < num_theta; j++ )
	    {
	    	 sintheta   = get_theta( num_theta, 0, j );
	    	 wtheta     = get_theta( num_theta, 1, j );
	    	 (*mu)[k]  	= cosphi*sintheta;
	    	 (*eta)[k]  = sinphi*sintheta;
	    	 d_wt[k]    = wphi*wtheta;
	    	 ++k;
	    }
	}
	Ensure ( k == num_ang / 4 ); // k must be the num of angles per quad

	// Octant directions
	double s[8][3] = { {1,1,1},{-1,1,1},{1,-1,1},{-1,-1,1},
			           {1,1,-1},{-1,1,-1},{1,-1,-1},{-1,-1,-1} };

	// Fill the next three quadrants
	for ( int o = 1; o < num_oct; o++ )
	{
		for ( int i = 0; i < num_ang/4; i++ )
		{

			(*mu)[k]  	= s[o][MU]*(*mu)[i] ;
			(*eta)[k] 	= s[o][ETA]*(*eta)[i] ;
			d_wt[k] 	= d_wt[i];
			++k;
		}
	}
	Ensure ( k == num_ang ); // k must now be the total number of angles

	// Normalize the quadrature set
	double wsum = 0.0;
	for (size_t n = 0; n <= num_ang - 1; ++n)
		wsum = wsum + d_wt[n];

	for (size_t n = 0; n <= num_ang - 1; ++n)
		d_wt[n] = d_wt[n] * (norm / wsum);

#ifdef DEBUG
	display();
#endif

	Ensure (d_angles[MU].size() == num_ang);
	Ensure (d_angles[ETA].size() == num_ang);

}

//---------------------------------------------------------------------------//
// PUBLIC FUNCTIONS
//---------------------------------------------------------------------------//
/*!
 * \brief Display the quadrature.
 */
void Quadruple_Range_2D::display() const
{

	using std::cout;
	using std::endl;
	using std::printf;
	Require ( d_angles.size()==2 );
	cout << endl << "The Quadrature directions and weights are:"
			<< endl << endl;
	cout << "   m         mu              eta             wt       " << endl;
	cout << "  ---   --------------  --------------  --------------" << endl;
	double sum_wt = 0.0;
	for ( int ix = 0; ix < d_angles[MU].size(); ++ix )
	{
		printf ("%4i    %13.10f   %13.10f   %13.10f   \n",
				ix, d_angles[MU][ix], d_angles[ETA][ix], d_wt[ix] );
		sum_wt += d_wt[ix];
	}
	cout << endl << "  The sum of the weights is " << sum_wt << endl;
	cout << endl;

}

/*
 * \brief return cos(phi) or weight
 *
 * \param N  number of azimuthal angles per quadrant
 * \param i  0 for angle and 1 for weight
 * \param j  index of requested point
 *
 */
double Quadruple_Range_2D::get_phi(int N, int i, int j)
{
	double a[30];
	double b[30];
	if( N == 2 )
	{
		a[0] = 0.328861319306339;
		a[1] = 0.9443782254288237;
		b[0] = 0.7853981633974483;
		b[1] = 0.7853981633974483;
	}
	else if( N == 4 )
	{
		a[0] = 0.1119194178021025;
		a[1] = 0.4990088153604911;
		a[2] = 0.8665969087139183;
		a[3] = 0.9937172857099943;
		b[0] = 0.2790495092691071;
		b[1] = 0.5063486541283412;
		b[2] = 0.5063486541283412;
		b[3] = 0.2790495092691071;
	}
	else if( N == 6 )
	{
		a[0] = 0.5498238339818787E-1;
		a[1] = 0.2686450605979331;
		a[2] = 0.5658713393447228;
		a[3] = 0.8244935580756284;
		a[4] = 0.9632392389309796;
		a[5] = 0.9984873246645922;
		b[0] = 0.1389912271015858;
		b[1] = 0.2849847975948065;
		b[2] = 0.361422138701056;
		b[3] = 0.361422138701056;
		b[4] = 0.2849847975948065;
		b[5] = 0.1389912271015858;
	}
	else if( N == 8 )
	{
		a[0] = 0.3248697884632087E-1;
		a[1] = 0.1639512569400686;
		a[2] = 0.370037938547133;
		a[3] = 0.6006186663088213;
		a[4] = 0.7995356262740346;
		a[5] = 0.9290166435730784;
		a[6] = 0.9864684411311756;
		a[7] = 0.9994721587945503;
		b[0] = 0.8261058789185801E-1;
		b[1] = 0.1780108196819361;
		b[2] = 0.2453650131486301;
		b[3] = 0.2794117426750241;
		b[4] = 0.2794117426750241;
		b[5] = 0.2453650131486301;
		b[6] = 0.1780108196819361;
		b[7] = 0.8261058789185801E-1;
	}
	else if( N == 10 )
	{
		a[0] = 0.2140636649836879E-1;
		a[1] = 0.1096782080789725;
		a[2] = 0.2554062770706209;
		a[3] = 0.4359965790249476;
		a[4] = 0.6217587837057514;
		a[5] = 0.7832087939271013;
		a[6] = 0.8999483224488741;
		a[7] = 0.9668338190366146;
		a[8] = 0.9939671476827519;
		a[9] = 0.9997708574835223;
		b[0] = 0.545988342771959E-1;
		b[1] = 0.1206689353273136;
		b[2] = 0.1734506495359569;
		b[3] = 0.2093517367248734;
		b[4] = 0.2273280075321085;
		b[5] = 0.2273280075321085;
		b[6] = 0.2093517367248734;
		b[7] = 0.1734506495359569;
		b[8] = 0.1206689353273136;
		b[9] = 0.545988342771959E-1;
	}
	else if( N == 12 )
	{
		a[0] = 0.1515427633810212E-1;
		a[1] = 0.7829331784823644E-1;
		a[2] = 0.1853856824160965;
		a[3] = 0.3250513900479526;
		a[4] = 0.481334101720309;
		a[5] = 0.6359378348511013;
		a[6] = 0.7717402867577236;
		a[7] = 0.8765372111445716;
		a[8] = 0.945696353925452;
		a[9] = 0.982665837787759;
		a[10] = 0.99693036686637;
		a[11] = 0.9998851673610667;
		b[0] = 0.3871995573990722E-1;
		b[1] = 0.8683898824818842E-1;
		b[2] = 0.1279220161079748;
		b[3] = 0.1595962888003237;
		b[4] = 0.1808538704864877;
		b[5] = 0.1914670440145664;
		b[6] = 0.1914670440145664;
		b[7] = 0.1808538704864877;
		b[8] = 0.1595962888003237;
		b[9] = 0.1279220161079748;
		b[10] = 0.8683898824818842E-1;
		b[11] = 0.3871995573990722E-1;
	}
	else if( N == 14 )
	{
		a[0] = 0.1128664599103795E-1;
		a[1] = 0.5860802825933166E-1;
		a[2] = 0.1401592977980301;
		a[3] = 0.2496535423349043;
		a[4] = 0.3780145294131805;
		a[5] = 0.5141044083621162;
		a[6] = 0.6460952816893911;
		a[7] = 0.7632567634673841;
		a[8] = 0.8577276125336286;
		a[9] = 0.92579966275244;
		a[10] = 0.9683352254253865;
		a[11] = 0.9901289669738802;
		a[12] = 0.9982810721553091;
		a[13] = 0.9999363037825324;
		b[0] = 0.2886965398857025E-1;
		b[1] = 0.6534816595720066E-1;
		b[2] = 0.9778151915986106E-1;
		b[3] = 0.124592860220795;
		b[4] = 0.1449498625863828;
		b[5] = 0.1585386998315389;
		b[6] = 0.1653174016530996;
		b[7] = 0.1653174016530996;
		b[8] = 0.1585386998315389;
		b[9] = 0.1449498625863828;
		b[10] = 0.124592860220795;
		b[11] = 0.9778151915986106E-1;
		b[12] = 0.6534816595720066E-1;
		b[13] = 0.2886965398857025E-1;
	}
	else if( N == 16 )
	{
		a[0] = 0.8729657082861998E-2;
		a[1] = 0.4548138640446708E-1;
		a[2] = 0.1094657179495637;
		a[3] = 0.1969448911359501;
		a[4] = 0.3024482742232713;
		a[5] = 0.4190211144820449;
		a[6] = 0.5387770134090762;
		a[7] = 0.6537245138927594;
		a[8] = 0.7567326211652801;
		a[9] = 0.8424484137453142;
		a[10] = 0.9079764895735049;
		a[11] = 0.9531657995434818;
		a[12] = 0.9804145602016775;
		a[13] = 0.993990571682542;
		a[14] = 0.9989651863256935;
		a[15] = 0.9999618958176435;
		b[0] = 0.2234563433970334E-1;
		b[1] = 0.5089471201185089E-1;
		b[2] = 0.7696546183192337E-1;
		b[3] = 0.9948832189094168E-1;
		b[4] = 0.1178163725528086;
		b[5] = 0.1316318495967909;
		b[6] = 0.1408332013024988;
		b[7] = 0.1454226098709307;
		b[8] = 0.1454226098709307;
		b[9] = 0.1408332013024988;
		b[10] = 0.1316318495967909;
		b[11] = 0.1178163725528086;
		b[12] = 0.9948832189094168E-1;
		b[13] = 0.7696546183192337E-1;
		b[14] = 0.5089471201185089E-1;
		b[15] = 0.2234563433970334E-1;
	}
	else if( N == 18 )
	{
		a[0] = 0.6952010366745025E-2;
		a[1] = 0.3630329822321955E-1;
		a[2] = 0.877584299640462E-1;
		a[3] = 0.1589616059026264;
		a[4] = 0.2464422434563339;
		a[5] = 0.3457022371739597;
		a[6] = 0.4514464196113868;
		a[7] = 0.5579698061333865;
		a[8] = 0.6596626820591664;
		a[9] = 0.7515618044435914;
		a[10] = 0.8298612507181373;
		a[11] = 0.8922982294166338;
		a[12] = 0.9383442668940432;
		a[13] = 0.9691574797937686;
		a[14] = 0.9872847653280476;
		a[15] = 0.996141786077788;
		a[16] = 0.9993408180086091;
		a[17] = 0.9999758344839443;
		b[0] = 0.1780448256851573E-1;
		b[1] = 0.4072837739441321E-1;
		b[2] = 0.620551874722286E-1;
		b[3] = 0.8104138647658412E-1;
		b[4] = 0.9719456818914889E-1;
		b[5] = 0.1102303994594047;
		b[6] = 0.1200213906455903;
		b[7] = 0.1265362696526889;
		b[8] = 0.1297861015388739;
		b[9] = 0.1297861015388739;
		b[10] = 0.1265362696526889;
		b[11] = 0.1200213906455903;
		b[12] = 0.1102303994594047;
		b[13] = 0.9719456818914889E-1;
		b[14] = 0.8104138647658412E-1;
		b[15] = 0.620551874722286E-1;
		b[16] = 0.4072837739441321E-1;
		b[17] = 0.1780448256851573E-1;
	}
	else if( N == 20 )
	{
		a[0] = 0.5666555435801935E-2;
		a[1] = 0.29639792746951E-1;
		a[2] = 0.7187424763201411E-1;
		a[3] = 0.1308135974586322;
		a[4] = 0.2041579376098809;
		a[5] = 0.2888941786344192;
		a[6] = 0.3813947484262091;
		a[7] = 0.4776033421571356;
		a[8] = 0.5732994237606388;
		a[9] = 0.6644147983089062;
		a[10] = 0.7473640182589308;
		a[11] = 0.8193459408062747;
		a[12] = 0.8785755787411428;
		a[13] = 0.9244122705118689;
		a[14] = 0.9573610361567596;
		a[15] = 0.9789379635660678;
		a[16] = 0.9914069813754243;
		a[17] = 0.997413701794462;
		a[18] = 0.9995606448264746;
		a[19] = 0.9999839449458641;
		b[0] = 0.145177950931735E-1;
		b[1] = 0.3331533618587398E-1;
		b[2] = 0.5104054858395153E-1;
		b[3] = 0.67163625891037E-1;
		b[4] = 0.8131097620559026E-1;
		b[5] = 0.9324211470049753E-1;
		b[6] = 0.1028266819400935;
		b[7] = 0.1100123318316945;
		b[8] = 0.1147920152181985;
		b[9] = 0.117176737747338;
		b[10] = 0.117176737747338;
		b[11] = 0.1147920152181985;
		b[12] = 0.1100123318316945;
		b[13] = 0.1028266819400935;
		b[14] = 0.9324211470049753E-1;
		b[15] = 0.8131097620559026E-1;
		b[16] = 0.67163625891037E-1;
		b[17] = 0.5104054858395153E-1;
		b[18] = 0.3331533618587398E-1;
		b[19] = 0.145177950931735E-1;
	}
	else if( N == 22 )
	{
		a[0] = 0.4707114717103505E-2;
		a[1] = 0.2465173376735076E-1;
		a[2] = 0.5991670409568439E-1;
		a[3] = 0.1094334822143285;
		a[4] = 0.1716206775350476;
		a[5] = 0.2443940481368066;
		a[6] = 0.3252092377487526;
		a[7] = 0.4111529714803756;
		a[8] = 0.4990835442262655;
		a[9] = 0.5858114200649882;
		a[10] = 0.668303327234281;
		a[11] = 0.7438888779969691;
		a[12] = 0.8104473950364958;
		a[13] = 0.8665538736181088;
		a[14] = 0.9115663629395599;
		a[15] = 0.945642084344217;
		a[16] = 0.9696760021962513;
		a[17] = 0.9851631047914915;
		a[18] = 0.9939941211951136;
		a[19] = 0.998203380364097;
		a[20] = 0.9996960998334762;
		a[21] = 0.9999889214741531;
		b[0] = 0.1206308084642689E-1;
		b[1] = 0.2774821193462954E-1;
		b[2] = 0.4268848001641819E-1;
		b[3] = 0.5649746567994538E-1;
		b[4] = 0.6889030401909524E-1;
		b[5] = 0.7966885558279603E-1;
		b[6] = 0.8871167329713092E-1;
		b[7] = 0.9595714699102763E-1;
		b[8] = 0.1013841131358227;
		b[9] = 0.1049937517313658;
		b[10] = 0.10679508016279;
		b[11] = 0.10679508016279;
		b[12] = 0.1049937517313658;
		b[13] = 0.1013841131358227;
		b[14] = 0.9595714699102763E-1;
		b[15] = 0.8871167329713092E-1;
		b[16] = 0.7966885558279603E-1;
		b[17] = 0.6889030401909524E-1;
		b[18] = 0.5649746567994538E-1;
		b[19] = 0.4268848001641819E-1;
		b[20] = 0.2774821193462954E-1;
		b[21] = 0.1206308084642689E-1;
	}
	else if( N == 24 )
	{
		a[0] = 0.3972127715469235E-2;
		a[1] = 0.208221968421092E-1;
		a[2] = 0.5069784852121022E-1;
		a[3] = 0.9284156235089204E-1;
		a[4] = 0.1461314233398479;
		a[5] = 0.2090839313632731;
		a[6] = 0.2798722573569018;
		a[7] = 0.3563718475936031;
		a[8] = 0.4362363542057792;
		a[9] = 0.5170013989086162;
		a[10] = 0.5962091754811982;
		a[11] = 0.671543718038901;
		a[12] = 0.740964934907509;
		a[13] = 0.802829134419043;
		a[14] = 0.855984552154146;
		a[15] = 0.8998321194918805;
		a[16] = 0.9343442118634447;
		a[17] = 0.9600372490491982;
		a[18] = 0.9778976989673706;
		a[19] = 0.9892651854344568;
		a[20] = 0.9956808948153246;
		a[21] = 0.9987140372275341;
		a[22] = 0.9997831945570342;
		a[23] = 0.9999921110695884;
		b[0] = 0.1018170331312922E-1;
		b[1] = 0.2346349353303599E-1;
		b[2] = 0.3621292193062666E-1;
		b[3] = 0.4814205277021821E-1;
		b[4] = 0.5903176230890425E-1;
		b[5] = 0.6872060753743357E-1;
		b[6] = 0.7710072080607362E-1;
		b[7] = 0.8410905852056142E-1;
		b[8] = 0.8971603507554962E-1;
		b[9] = 0.9391384633073917E-1;
		b[10] = 0.9670612366917621E-1;
		b[11] = 0.9809983760200038E-1;
		b[12] = 0.9809983760200038E-1;
		b[13] = 0.9670612366917621E-1;
		b[14] = 0.9391384633073917E-1;
		b[15] = 0.8971603507554962E-1;
		b[16] = 0.8410905852056142E-1;
		b[17] = 0.7710072080607362E-1;
		b[18] = 0.6872060753743357E-1;
		b[19] = 0.5903176230890425E-1;
		b[20] = 0.4814205277021821E-1;
		b[21] = 0.3621292193062666E-1;
		b[22] = 0.2346349353303599E-1;
		b[23] = 0.1018170331312922E-1;
	}
	else if( N == 26 )
	{
		a[0] = 0.3396696521564986E-2;
		a[1] = 0.1781892050929334E-1;
		a[2] = 0.4344484549656114E-1;
		a[3] = 0.7972274758288012E-1;
		a[4] = 0.1258351250930473;
		a[5] = 0.180699198021927;
		a[6] = 0.2429742598458005;
		a[7] = 0.3110847255339362;
		a[8] = 0.3832618562348258;
		a[9] = 0.4576040998675997;
		a[10] = 0.5321533124724104;
		a[11] = 0.6049820056372407;
		a[12] = 0.6742853437025943;
		a[13] = 0.7384709034673433;
		a[14] = 0.7962391430061333;
		a[15] = 0.8466480095202735;
		a[16] = 0.8891560536735741;
		a[17] = 0.9236397293076104;
		a[18] = 0.950382182881958;
		a[19] = 0.9700327360725438;
		a[20] = 0.9835384079100482;
		a[21] = 0.9920511686867856;
		a[22] = 0.9968170762571418;
		a[23] = 0.9990558269685333;
		a[24] = 0.9998412304320539;
		a[25] = 0.9999942312097307;
		b[0] = 0.870818485865091E-2;
		b[1] = 0.2009672209559701E-1;
		b[2] = 0.3109543899569998E-1;
		b[3] = 0.4148557359431433E-1;
		b[4] = 0.510968162469414E-1;
		b[5] = 0.5979803209497304E-1;
		b[6] = 0.6749583467970308E-1;
		b[7] = 0.7413013546424966E-1;
		b[8] = 0.7966752788841578E-1;
		b[9] = 0.8409386224353523E-1;
		b[10] = 0.8740711539780842E-1;
		b[11] = 0.8961127788931235E-1;
		b[12] = 0.9071164194824712E-1;
		b[13] = 0.9071164194824712E-1;
		b[14] = 0.8961127788931235E-1;
		b[15] = 0.8740711539780842E-1;
		b[16] = 0.8409386224353523E-1;
		b[17] = 0.7966752788841578E-1;
		b[18] = 0.7413013546424966E-1;
		b[19] = 0.6749583467970308E-1;
		b[20] = 0.5979803209497304E-1;
		b[21] = 0.510968162469414E-1;
		b[22] = 0.4148557359431433E-1;
		b[23] = 0.3109543899569998E-1;
		b[24] = 0.2009672209559701E-1;
		b[25] = 0.870818485865091E-2;
	}
	else if( N == 28 )
	{
		a[0] = 0.2937783068161103E-2;
		a[1] = 0.1542057328035316E-1;
		a[2] = 0.3763822344160284E-1;
		a[3] = 0.6917953054276372E-1;
		a[4] = 0.1094349422308658;
		a[5] = 0.1575965218033132;
		a[6] = 0.2126604305088603;
		a[7] = 0.2734385687346029;
		a[8] = 0.3385818199466888;
		a[9] = 0.4066156050266695;
		a[10] = 0.4759868689664975;
		a[11] = 0.5451202571507631;
		a[12] = 0.6124802047142074;
		a[13] = 0.6766350014909858;
		a[14] = 0.736318595960535;
		a[15] = 0.7904859257654386;
		a[16] = 0.8383578622783268;
		a[17] = 0.8794523867563669;
		a[18] = 0.9135993376468678;
		a[19] = 0.9409369538930799;
		a[20] = 0.9618894682490145;
		a[21] = 0.9771261644720124;
		a[22] = 0.9875035880013286;
		a[23] = 0.9939939604539492;
		a[24] = 0.9976042264113975;
		a[25] = 0.9992914310330896;
		a[26] = 0.999881095890759;
		a[27] = 0.9999956847060113;
		b[0] = 0.7532680050608221E-2;
		b[1] = 0.1740393870836662E-1;
		b[2] = 0.2698374387900945E-1;
		b[3] = 0.3610294286409177E-1;
		b[4] = 0.4462772409676748E-1;
		b[5] = 0.5245143079823476E-1;
		b[6] = 0.5949430183433835E-1;
		b[7] = 0.6570130151997985E-1;
		b[8] = 0.7103829899589742E-1;
		b[9] = 0.7548738262548303E-1;
		b[10] = 0.7904202969199018E-1;
		b[11] = 0.8170266224145589E-1;
		b[12] = 0.834729185573813E-1;
		b[13] = 0.8435680753384398E-1;
		b[14] = 0.8435680753384398E-1;
		b[15] = 0.834729185573813E-1;
		b[16] = 0.8170266224145589E-1;
		b[17] = 0.7904202969199018E-1;
		b[18] = 0.7548738262548303E-1;
		b[19] = 0.7103829899589742E-1;
		b[20] = 0.6570130151997985E-1;
		b[21] = 0.5949430183433835E-1;
		b[22] = 0.5245143079823476E-1;
		b[23] = 0.4462772409676748E-1;
		b[24] = 0.3610294286409177E-1;
		b[25] = 0.2698374387900945E-1;
		b[26] = 0.1740393870836662E-1;
		b[27] = 0.7532680050608221E-2;
	}
	else if( N == 30 )
	{
		a[0] = 0.2565935816040464E-2;
		a[1] = 0.1347516554782726E-1;
		a[2] = 0.3291878087349994E-1;
		a[3] = 0.6058415580193666E-1;
		a[4] = 0.960076213078193E-1;
		a[5] = 0.1385744973107413;
		a[6] = 0.1875194470293309;
		a[7] = 0.2419322530851564;
		a[8] = 0.300770887782924;
		a[9] = 0.3628827241348177;
		a[10] = 0.4270337881301431;
		a[11] = 0.4919450854898168;
		a[12] = 0.5563343213452209;
		a[13] = 0.6189608105599492;
		a[14] = 0.6786710589083211;
		a[15] = 0.7344423692845193;
		a[16] = 0.7854218707108751;
		a[17] = 0.8309585566641411;
		a[18] = 0.870626230286233;
		a[19] = 0.9042356682830091;
		a[20] = 0.9318348182615274;
		a[21] = 0.9536964260509063;
		a[22] = 0.9702931438061076;
		a[23] = 0.9822608904897995;
		a[24] = 0.9903520125162948;
		a[25] = 0.9953805988921094;
		a[26] = 0.9981630929190713;
		a[27] = 0.9994580300671972;
		a[28] = 0.9999092058349391;
		a[29] = 0.9999967079812753;
		b[0] = 0.6579962619041378E-2;
		b[1] = 0.1521696776845086E-1;
		b[2] = 0.2363209088631991E-1;
		b[3] = 0.3169244521291083E-1;
		b[4] = 0.3929176131897425E-1;
		b[5] = 0.4634296332807871E-1;
		b[6] = 0.5277837509680457E-1;
		b[7] = 0.5854875269926412E-1;
		b[8] = 0.6362110627094798E-1;
		b[9] = 0.6797575914341677E-1;
		b[10] = 0.7160310605662001E-1;
		b[11] = 0.745004453275371E-1;
		b[12] = 0.7666914603657626E-1;
		b[13] = 0.781123059148995E-1;
		b[14] = 0.7883297571760606E-1;
		b[15] = 0.7883297571760606E-1;
		b[16] = 0.781123059148995E-1;
		b[17] = 0.7666914603657626E-1;
		b[18] = 0.745004453275371E-1;
		b[19] = 0.7160310605662001E-1;
		b[20] = 0.6797575914341677E-1;
		b[21] = 0.6362110627094798E-1;
		b[22] = 0.5854875269926412E-1;
		b[23] = 0.5277837509680457E-1;
		b[24] = 0.4634296332807871E-1;
		b[25] = 0.3929176131897425E-1;
		b[26] = 0.3169244521291083E-1;
		b[27] = 0.2363209088631991E-1;
		b[28] = 0.1521696776845086E-1;
		b[29] = 0.6579962619041378E-2;
	}
	else
	{
		throw util::assertion(
				"Unsupported value of N for phi.");
	}
	if ( i == 0 )
		return a[j]; // return cos(phi_j)
	else
		return b[j]; // return w_j

}

/*
 * \brief return sin(theta) or weight
 *
 * \param N  number of polar angles per quadrant
 * \param i  0 for angle and 1 for weight
 * \param j  index of requested point
 *
 */
double Quadruple_Range_2D::get_theta(int N, int i, int j)
{
	double a[15];
	double b[15];
	if( N == 1 )
	{
		a[0] = 0.7853981633974483;
		b[0] = 0.2E1;
	}
	else if( N == 2 )
	{
		a[0] = 0.3993735638823956;
		a[1] = 0.9144479484659753;
		b[0] = 0.5010918381152245;
		b[1] = 0.1498908161884776E1;
	}
	else if( N == 3 )
	{
		a[0] = 0.2311522961278847;
		a[1] = 0.6399711962340753;
		a[2] = 0.9544967673433067;
		b[0] = 0.1706016479452537;
		b[1] = 0.6829125158656247;
		b[2] = 0.1146485836189122E1;
	}
	else if( N == 4 )
	{
		a[0] = 0.149249146573134;
		a[1] = 0.4452174134055993;
		a[2] = 0.7657198751837065;
		a[3] = 0.971840383889606;
		b[0] = 0.7210020864182995E-1;
		b[1] = 0.3302019567476854;
		b[2] = 0.6776788331337992;
		b[3] = 0.9200190014766855;
	}
	else if( N == 5 )
	{
		a[0] = 0.1039469054750266;
		a[1] = 0.3218997600306228;
		a[2] = 0.5925641316724047;
		a[3] = 0.8367574981862629;
		a[4] = 0.9808766137326722;
		b[0] = 0.3530233934872431E-1;
		b[1] = 0.1745260006417964;
		b[2] = 0.3987720643885643;
		b[3] = 0.6254719597163901;
		b[4] = 0.7659276359045249;
	}
	else if( N == 6 )
	{
		a[0] = 0.764222577018977E-1;
		a[1] = 0.2418380224706412;
		a[2] = 0.4624149310771195;
		a[3] = 0.69176378901922;
		a[4] = 0.8801887364309769;
		a[5] = 0.9861717557161557;
		b[0] = 0.1920436605371306E-1;
		b[1] = 0.9964572033500919E-1;
		b[2] = 0.243394408480036;
		b[3] = 0.415494935541075;
		b[4] = 0.5670894260317977;
		b[5] = 0.655171143558369;
	}
	else if( N == 7 )
	{
		a[0] = 0.5849788427341292E-1;
		a[1] = 0.1876744099204198;
		a[2] = 0.3673808738342405;
		a[3] = 0.5691413298688445;
		a[4] = 0.76007398923149;
		a[5] = 0.9084990189032058;
		a[6] = 0.9895378462166649;
		b[0] = 0.1130288664498013E-1;
		b[1] = 0.6059668882540333E-1;
		b[2] = 0.1548577023350282;
		b[3] = 0.2795356837349958;
		b[4] = 0.408359039083128;
		b[5] = 0.5133389971363539;
		b[6] = 0.5720090022401106;
	}
	else if( N == 8 )
	{
		a[0] = 0.4619235546732977E-1;
		a[1] = 0.1495785939929699;
		a[2] = 0.2974141685414306;
		a[3] = 0.4714033535759624;
		a[4] = 0.6494576043472951;
		a[5] = 0.8085613275367436;
		a[6] = 0.9279142834477891;
		a[7] = 0.9918094732391086;
		b[0] = 0.7070596107672168E-2;
		b[1] = 0.3879562306402765E-1;
		b[2] = 0.1024016241431204;
		b[3] = 0.1923042113259944;
		b[4] = 0.2943821042432288;
		b[5] = 0.3912388348106103;
		b[6] = 0.4664219853334316;
		b[7] = 0.5073850209719148;
	}
	else if( N == 9 )
	{
		a[0] = 0.373871754313704E-1;
		a[1] = 0.1218673184482112;
		a[2] = 0.2449791366335987;
		a[3] = 0.3944837657444153;
		a[4] = 0.5551788717532553;
		a[5] = 0.7104369578733888;
		a[6] = 0.8440015993098882;
		a[7] = 0.9417816602325884;
		a[8] = 0.9934143180700493;
		b[0] = 0.4643065076049817E-2;
		b[1] = 0.259137055013675E-1;
		b[2] = 0.7006244894417409E-1;
		b[3] = 0.1354951021927222;
		b[4] = 0.2146387592743282;
		b[5] = 0.2967844729635604;
		b[6] = 0.3706102754254009;
		b[7] = 0.4260782931411163;
		b[8] = 0.4557738774812805;
	}
	else if( N == 10 )
	{
		a[0] = 0.3087344674320983E-1;
		a[1] = 0.1011256691908164;
		a[2] = 0.2049120334874169;
		a[3] = 0.3337559301018391;
		a[4] = 0.47692469846783;
		a[5] = 0.6222949084152428;
		a[6] = 0.7574090067646704;
		a[7] = 0.8705939190886369;
		a[8] = 0.9520192454053205;
		a[9] = 0.9945898978711666;
		b[0] = 0.317192548389715E-2;
		b[1] = 0.1793256207590737E-1;
		b[2] = 0.4938267089469711E-1;
		b[3] = 0.9768026105928111E-1;
		b[4] = 0.1588277298633008;
		b[5] = 0.2262349688595732;
		b[6] = 0.29226033258076;
		b[7] = 0.3494452142568355;
		b[8] = 0.39143378370197;
		b[9] = 0.4136305512237778;
	}
	else if( N == 11 )
	{
		a[0] = 0.259214348332132E-1;
		a[1] = 0.8521935262972822E-1;
		a[2] = 0.1737212859610505;
		a[3] = 0.2853788554402115;
		a[4] = 0.4124302930568668;
		a[5] = 0.545908389315292;
		a[6] = 0.6762815928489082;
		a[7] = 0.7941615804130648;
		a[8] = 0.8910104077090977;
		a[9] = 0.9597863457139877;
		a[10] = 0.9954766828497256;
		b[0] = 0.2239168223242646E-2;
		b[1] = 0.1278609020279452E-1;
		b[2] = 0.3571984407731606E-1;
		b[3] = 0.7192019606102708E-1;
		b[4] = 0.1193656258066426;
		b[5] = 0.1740015513470782;
		b[6] = 0.2306960076947862;
		b[7] = 0.2840568417450768;
		b[8] = 0.3290675185193199;
		b[9] = 0.3615644996448663;
		b[10] = 0.3785826566778497;
	}
	else if( N == 12 )
	{
		a[0] = 0.2206987544161146E-1;
		a[1] = 0.7276548794207686E-1;
		a[2] = 0.1490240471276492;
		a[3] = 0.2464178213464814;
		a[4] = 0.3592090638180165;
		a[5] = 0.4806598424300749;
		a[6] = 0.6034329607303313;
		a[7] = 0.7200483966382738;
		a[8] = 0.8233578535894617;
		a[9] = 0.9070015920820176;
		a[10] = 0.9658155920819829;
		a[11] = 0.9961620595917704;
		b[0] = 0.1625008098812313E-2;
		b[1] = 0.9352495000167818E-2;
		b[2] = 0.2642849768723387E-1;
		b[3] = 0.5397608657841666E-1;
		b[4] = 0.9107347191808735E-1;
		b[5] = 0.1352375815741906;
		b[6] = 0.1830231517042723;
		b[7] = 0.2305709517590055;
		b[8] = 0.2740530323427262;
		b[9] = 0.3100184495077959;
		b[10] = 0.3356558712699347;
		b[11] = 0.3489854025593568;
	}
	else if( N == 13 )
	{
		a[0] = 0.190157862781181E-1;
		a[1] = 0.6283891176988711E-1;
		a[2] = 0.1291669489733978;
		a[3] = 0.2146859046048244;
		a[4] = 0.3150698520346282;
		a[5] = 0.425176468911348;
		a[6] = 0.5393050111837124;
		a[7] = 0.6514970179728065;
		a[8] = 0.7558583950754834;
		a[9] = 0.846881415042289;
		a[10] = 0.9197467199541607;
		a[11] = 0.9705877364229169;
		a[12] = 0.9967027009889732;
		b[0] = 0.1207468708743215E-2;
		b[1] = 0.6993516363714714E-2;
		b[2] = 0.1994654540921309E-1;
		b[3] = 0.4121379610227348E-1;
		b[4] = 0.7048445872346682E-1;
		b[5] = 0.1062573893438944;
		b[6] = 0.1462194120104745;
		b[7] = 0.1876167252194712;
		b[8] = 0.2275704275088438;
		b[9] = 0.2633290245190361;
		b[10] = 0.2924659823505597;
		b[11] = 0.3130314389266094;
		b[12] = 0.3236638148136996;
	}
	else if( N == 14 )
	{
		a[0] = 0.1655357842785945E-1;
		a[1] = 0.5480298212338539E-1;
		a[2] = 0.1129813228893869;
		a[3] = 0.1885598024182952;
		a[4] = 0.2782189669515478;
		a[5] = 0.3779758905459106;
		a[6] = 0.483354614232666;
		a[7] = 0.5895884420495059;
		a[8] = 0.6918416750813907;
		a[9] = 0.7854377100908073;
		a[10] = 0.8660809241804746;
		a[11] = 0.9300608506655752;
		a[12] = 0.9744285442429917;
		a[13] = 0.9971366642642528;
		b[0] = 0.9156938974039617E-3;
		b[1] = 0.5331012020952745E-2;
		b[2] = 0.1532101378791769E-1;
		b[3] = 0.3196191025679207E-1;
		b[4] = 0.5527736741151233E-1;
		b[5] = 0.8438398073696996E-1;
		b[6] = 0.1177312073753468;
		b[7] = 0.1533547657512768;
		b[8] = 0.1891023104848334;
		b[9] = 0.2228197138483984;
		b[10] = 0.2525002447442476;
		b[11] = 0.2764024505279778;
		b[12] = 0.293141514050202;
		b[13] = 0.3017568151061684;
	}
	else if( N == 15 )
	{
		a[0] = 0.1453980971214771E-1;
		a[1] = 0.4820858869899936E-1;
		a[2] = 0.9962575637202815E-1;
		a[3] = 0.1668286525539518;
		a[4] = 0.247229707903158;
		a[5] = 0.3377013939031784;
		a[6] = 0.4346915962513905;
		a[7] = 0.5343624719763381;
		a[8] = 0.6327454393360694;
		a[9] = 0.7259042142853639;
		a[10] = 0.8100978365440864;
		a[11] = 0.8819361061028814;
		a[12] = 0.9385205574600149;
		a[13] = 0.9775649540709922;
		a[14] = 0.9974902738450473;
		b[0] = 0.7068830110205541E-3;
		b[1] = 0.4132909575263126E-2;
		b[2] = 0.1195308097311956E-1;
		b[3] = 0.2513687784360474E-1;
		b[4] = 0.4388448158131325E-1;
		b[5] = 0.6770273244117886E-1;
		b[6] = 0.9555753069484589E-1;
		b[7] = 0.1260476011692554;
		b[8] = 0.1575670106423901;
		b[9] = 0.1884437898269261;
		b[10] = 0.2170534569913444;
		b[11] = 0.2419105398978403;
		b[12] = 0.2617415950418412;
		b[13] = 0.2755422071013228;
		b[14] = 0.2826193032087338;
	}
	else
	{
		throw util::assertion(
				"Unsupported value of N for theta.");
	}
	if ( i==0 )
		return a[j]; // return sin(theta_j)
	else
		return b[j]; // return w_j

}


} // end namespace slabtran

//---------------------------------------------------------------------------//
//                 end of Quadruple_Range_2D.cc
//---------------------------------------------------------------------------//
